Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for effects plots in ggplotjk
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(DHARMa) # for assessing dispersion etc
library(glmmTMB) # for glmmTMB
library(performance) # for diagnostic plots
library(see) # for diagnostic plots
library(patchwork) # for multiple plots

Scenario

fanworms

In an attempt to understand the effects on marine animals of short-term exposure to toxic substances, such as might occur following a spill, or a major increase in storm water flows, a it was decided to examine the toxicant in question, Copper, as part of a field experiment in Hong Kong. The experiment consisted of small sources of Cu (small, hemispherical plaster blocks, impregnated with copper), which released the metal into sea water over 4 or 5 days. The organism whose response to Cu was being measured was a small, polychaete worm, Hydroides, that attaches to hard surfaces in the sea, and is one of the first species to colonise any surface that is submerged. The biological questions focused on whether the timing of exposure to Cu affects the overall abundance of these worms. The time period of interest was the first or second week after a surface being available.

The experimental setup consisted of sheets of black perspex (settlement plates), which provided good surfaces for these worms. Each plate had a plaster block bolted to its centre, and the dissolving block would create a gradient of [Cu] across the plate. Over the two weeks of the experiment, a given plate would have plain plaster blocks (Control) or a block containing copper in the first week, followed by a plain block, or a plain block in the first week, followed by a dose of copper in the second week. After two weeks in the water, plates were removed and counted back in the laboratory. Without a clear idea of how sensitive these worms are to copper, an effect of the treatments might show up as an overall difference in the density of worms across a plate, or it could show up as a gradient in abundance across the plate, with a different gradient in different treatments. Therefore, on each plate, the density of worms (#/cm2) was recorded at each of four distances from the centre of the plate.

Format of copper.csv data file

COPPER PLATE DIST WORMS AREA COUNT
.. .. .. .. .. ..
COPPER Categorical listing of the copper treatment (control = no copper applied, week 2 = copper treatment applied in second week and week 1= copper treatment applied in first week) applied to whole plates. Factor A (between plot factor).
PLATE Substrate provided for polychaete worm colonisation on which copper treatment applied. These are the plots (Factor B). Numbers in this column represent numerical labels given to each plate.
DIST Categorical listing for the four concentric distances from the center of the plate (source of copper treatment) with 1 being the closest and 4 the furthest. Factor C (within plot factor)
WORMS Density (#/cm2) of worms measured. Response variable.

Read in the data

copper <- read_csv("../public/data/copper.csv", trim_ws = TRUE)
glimpse(copper)
## Rows: 60
## Columns: 4
## $ COPPER <chr> "control", "control", "control", "control", "control", "control…
## $ PLATE  <dbl> 200, 200, 200, 200, 39, 39, 39, 39, 34, 34, 34, 34, 36, 36, 36,…
## $ DIST   <dbl> 4, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 1, 4, 3, 2, 1, 4, …
## $ WORMS  <dbl> 11.50, 13.00, 13.50, 12.00, 17.75, 13.75, 12.75, 12.50, 11.50, …

Data preparation

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{Gamma}(\lambda_i, \phi)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of copper, distance and their interaction on the number of number of worms. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual plates. \(\phi\) represents the shape parameter of the Gamma distribution.

Fit the model

Model validation

Partial plots

Model Investigation / hypothesis testing

Further analysis

Summary figures

newdata <- copper.glmmTMB4a %>%
  emmeans(~ COPPER | DIST, type = "response") %>%
  summary(infer = TRUE) %>%
  as.data.frame()

head(newdata)
g1 <- ggplot(newdata, aes(y = response, x = DIST, fill = COPPER)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),
    shape = 21,
    position = position_dodge(width = 0.2)
  ) +
  theme_bw()
copper.comp <- copper.glmmTMB4a %>%
  emmeans(~ COPPER | DIST, type = "response") %>%
  pairs() %>%
  ## contrast(method='pairwise') %>%
  summary(infer = TRUE) %>%
  as.data.frame()
head(copper.comp)
g2 <- ggplot(copper.comp, aes(y = ratio, x = contrast, color = DIST)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6)
  ) +
  geom_hline(yintercept = 1) +
  scale_y_continuous(trans = scales::log2_trans(), limits = c(0.25, 500)) +
  coord_flip() +
  theme_bw()

g1 + g2

References